﻿using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;

namespace DynamicProgramming
{
    /// <summary>
    /// Dynamic programming example
    /// Needleman-Wunsch algorithm
    /// Global sequence alignment using a scoring matrix (BLOSUM62)
    /// 
    /// Author:  Kunuk Nykjaer
    /// </summary>
    
    public class SequenceAlignment
    {
        public static void Main(string[] args)
        {
            new SequenceAlignment().Run();
            PL("Press a key to exit.");
            Console.ReadKey();
        }

        const int Size = 24; // BLOSUM62 matrix
        readonly private List<string> linesMatrix;
        readonly List<LifeForm> lifeForms = new List<LifeForm>();
        readonly private int[,] matrix = new int[Size, Size];
        readonly Dictionary<string, int> lookup = new Dictionary<string, int>();
        const int DefaultMismatchPenalty = -4; // if not defined, use this
        static readonly string NL = Environment.NewLine;

        // trace back
        const string DONE = @"¤";
        const string DIAG = @"\";
        const string UP = @"|";
        const string LEFT = @"-";

        // print alignment
        const string GAP = @"-";

        public SequenceAlignment()
        {
            linesMatrix = ReadFileToList("BLOSUM62.txt");
        }


        // BLOSUM62.txt
        void ParseMatrixFile()
        {
            const int startLine = 7; // # lines before line 7
            for (int i = 0; i < linesMatrix.Count; i++)
            {
                if (i >= startLine)
                {
                    var line = linesMatrix.ElementAt(i);
                    String[] row = Regex.Split(line, @"\s+"); // white space
                    lookup.Add(row[0], i - startLine);

                    for (int j = 1; j < row.Length; j++)
                    {
                        var s = row[j];
                        if (string.IsNullOrEmpty(s))
                            continue;
                        matrix[i - startLine, j - 1] = Int32.Parse(s);
                    }
                }
            }
        }

        static void PrintMatrix<T>(T[,] A, int I, int J)
        {
            for (int i = 0; i < I; i++)
            {
                for (int j = 0; j < J; j++)
                {
                    var v = A[i, j];
                    P(v + " ");
                }
                PL();
            }
        }


        static string ReverseString(string s)
        {
            char[] arr = s.ToCharArray();
            Array.Reverse(arr);
            return new string(arr);
        }

        public void Run()
        {
            // put your data here
            lifeForms.Add(new LifeForm() { Name = "Sphinx", DNA = "KQRK" });
            lifeForms.Add(new LifeForm() { Name = "Bandersnatch", DNA = "KAK" });
            lifeForms.Add(new LifeForm() { Name = "Snark", DNA = "KQRIKAAKABK" });

            ParseMatrixFile();

            for (int i = 0; i < lifeForms.Count; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    if (i == j)
                        continue;

                    var l1 = lifeForms.ElementAt(i);
                    var l2 = lifeForms.ElementAt(j);
                    var sequenceAlign = SequenceAlign(l1.DNA, l2.DNA).ToString();
                    Console.WriteLine(string.Format("Matching species: {0} and {1}, DNA: {2} and {3}\n{4}"
                        , l1.Name, l2.Name, l1.DNA, l2.DNA, sequenceAlign));
                }
            }
        }


        /*  
             p = gap penalty      
             a = diff char penalty
         
             Sequence-Alignment(m, n, xs, ys, p, a) {
             for i = 0 to m
                 M[i, 0] = i*p
             for j = 0 to n
                 M[0, j] = j*p
             for i = 1 to m
                 for j = 1 to n
                     M[i, j] = max(  a(xi, yj) + M[i-1, j-1],
                         p + M[i-1, j],
                         p + M[i, j-1])
             return M[m, n]            
         */
        Sequence SequenceAlign(string xs, string ys)
        {
            const int p = -4; //gap penalty, knowledge by looking at matrix file
            int m = xs.Length;
            int n = ys.Length;

            // init the matrix
            var M = new int[m + 1, n + 1]; // dynamic programming buttom up memory table
            var T = new string[m + 1, n + 1]; // trace back

            for (int i = 0; i < m + 1; i++)
                M[i, 0] = i * p;
            for (int j = 0; j < n + 1; j++)
                M[0, j] = j * p;

            T[0, 0] = DONE;
            for (int i = 1; i < m + 1; i++)
                T[i, 0] = UP;
            for (int j = 1; j < n + 1; j++)
                T[0, j] = LEFT;

            // calc
            for (int i = 1; i < m + 1; i++)
            {
                for (int j = 1; j < n + 1; j++)
                {
                    var alpha = Alpha(xs.ElementAt(i - 1).ToString(), ys.ElementAt(j - 1).ToString());
                    var diag = alpha + M[i - 1, j - 1];
                    var up = p + M[i - 1, j];
                    var left = p + M[i, j - 1];
                    var max = Max(diag, up, left);
                    M[i, j] = max;

                    if (max == diag)
                        T[i, j] = DIAG;
                    else if (max == up)
                        T[i, j] = UP;
                    else
                        T[i, j] = LEFT;
                }
            }

            var traceBack = ParseTraceBack(T, m + 1, n + 1);

            var sb = new StringBuilder();
            string first, second;

            if (xs.Length != ys.Length)
            {
                string s;
                if (xs.Length > ys.Length)
                {
                    s = ys;
                    first = xs;
                }
                else
                {
                    s = xs;
                    first = ys;
                }


                int i = 0;
                foreach (var trace in traceBack)
                {
                    if (trace.ToString() == DIAG)
                        sb.Append(s.ElementAt(i++).ToString());
                    else
                        sb.Append(GAP);
                }

                second = sb.ToString();
            }
            else
            {
                first = xs;
                second = ys;
            }

            PL("\nScore table");
            PrintMatrix(M, m + 1, n + 1);
            PL("\nTraceBack");
            PrintMatrix(T, m + 1, n + 1);
            PL();

            var sequence = new Sequence() { Score = M[m, n], Path = traceBack, One = first, Two = second };
            return sequence;
        }


        static string ParseTraceBack(string[,] T, int I, int J)
        {
            var sb = new StringBuilder();
            int i = I - 1;
            int j = J - 1;
            string path = T[i, j];
            while (path != DONE)
            {
                sb.Append(path);
                if (path == DIAG)
                {
                    i--;
                    j--;
                }
                else if (path == UP)
                    i--;
                else if (path == LEFT)
                    j--;

                path = T[i, j];
            }
            return ReverseString(sb.ToString());
        }


        static int Max(int a, int b, int c)
        {
            if (a >= b && a >= c)
                return a;
            if (b >= a && b >= c)
                return b;
            return c;
        }

        int Alpha(string x, string y)
        {
            if (lookup.ContainsKey(x) && lookup.ContainsKey(y))
            {
                var i = lookup[x];
                var j = lookup[y];
                return matrix[i, j];
            }
            else if (x == y)
                return 1; // matrix file match * with *

            return DefaultMismatchPenalty; // default mismatch penalty
        }

        class LifeForm
        {
            public string Name { get; set; }
            public string DNA { get; set; }
        }

        class Sequence
        {
            public int Score { get; set; }
            public string Path { get; set; }
            public string One { get; set; }
            public string Two { get; set; }
            public new string ToString()
            {
                var s = string.Format("score = {0}{1}one = {2}{3}two = {4}\n\n", Score, NL, One, NL, Two);
                return s;
            }
        }

        public static List<string> ReadFileToList(string path)
        {
            var lines = new List<string>();
            StreamReader streamReader = null;
            try
            {
                streamReader = File.OpenText(path);
                string line = streamReader.ReadLine();
                while (line != null)
                {
                    lines.Add(line);
                    line = streamReader.ReadLine();
                }
            }
            catch (FileNotFoundException ex)
            {
                PL(string.Format("Error cannot find file {0}", path));
                PL(ex.Message + Environment.NewLine + ex.StackTrace);
            }
            catch (Exception ex)
            {
                PL(string.Format("Error reading file {0}", path));
                PL(ex.Message + Environment.NewLine + ex.StackTrace);
            }
            finally
            {
                if (streamReader != null)
                    streamReader.Close();
            }
            return lines;
        }

        public static void PL(object o) { Console.WriteLine(o); } //alias
        public static void PL() { Console.WriteLine(); } //alias
        public static void P(object o) { Console.Write(o); } //alias
    }
}
